library(MASS)
#library(Zelig)
library(gplots)
library(ggplot2)
library(Hmisc)
library(sciplot)
library(lme4)
library(arm)
library(calibrate)
library(memisc)

rm(list=ls())

### This is the setup for Figure 4


data <- read.csv(file="working_2014_07_11.csv",head=TRUE,sep=",")

n <- length(data$logimports96)
y <- data$logimports96
x <- data$dispute

st1hsnum <- data$st1hs_num
stnam1 <- data$stnam1
year <- data$year
countryyear<-data$st1_year_num

### Figure 4: Country specific dispute coefficients
# x|country: x is only the dispute dummy variables
# 1|countryyear, 1|st1hsnum: the d_it and a_is random effects
m2cyfe <- lmer(y ~ (1|countryyear) + (x|stnam1) + (1|st1hsnum) )
m2cyfe

#test <- lmer(y ~ x + (1|countryyear) + (x|stnam1) + (1|st1hsnum) )
#test

### Making the figure with country coefficients
b.hat.m2b <- ranef(m2cyfe)$stnam1[,2]
b.se.m2b <- se.ranef(m2cyfe )$stnam1[,2]
J <- 31
countryvec <- as.numeric(1:J)

# Make vector of country names
data.nomissname <- data[ which(data$stnam1 != ""), ]
country.name <- as.vector(data.nomissname$stnam1) ### Need to remove missing names
uniq.name <- sort(unique(country.name))

lowerb <- as.matrix(b.hat.m2b - 2*b.se.m2b, nrow=J, ncol=1)
upperb <- as.matrix(b.hat.m2b + 2*b.se.m2b, nrow=J, ncol=1)
b.hatb <- as.matrix(b.hat.m2b, nrow=J, ncol=1)
country.num.unorderb <- as.matrix(countryvec, nrow=J, ncol=1)
country.nam.unorderb <- as.matrix(uniq.name, nrow=J, ncol=1)

fig.matb <- as.data.frame(cbind(country.num.unorderb, b.hatb, lowerb, upperb))
fig.matb <- as.data.frame(cbind(fig.matb, country.nam.unorderb))

fig.matb <- fig.matb[order(fig.matb[,2]),]

country.num.orderb <- as.numeric(1:J)
country.num.orderb <- as.matrix(country.num.orderb, nrow=J, ncol=1)
fig.matb <- as.data.frame(cbind(fig.matb, country.num.orderb))


# Making the Figure
plot(fig.matb$country.num.order, fig.matb$V2, xlab= expression(paste("Country i, ranked by ", beta[i])), ylab= expression(beta[i]))
segments(fig.matb$country.num.order, fig.matb$V3, fig.matb$country.num.order, fig.matb$V4, lwd=.5, col="gray10")
abline(h=0)
text(fig.matb$country.num.order, fig.matb$V2, fig.matb$country.nam.unorder, cex=0.6, pos=3, offset = 0.48,  col="red") 

# Making the Figure wider
plot(fig.matb$country.num.order, fig.matb$V2, xlab= expression(paste("Country i, ranked by ", beta[i])), ylab= expression(beta[i]))
segments(fig.matb$country.num.order, fig.matb$V3, fig.matb$country.num.order, fig.matb$V4, lwd=.5, col="gray10")
abline(h=0)
text(fig.matb$country.num.order, fig.matb$V2, fig.matb$country.nam.unorder, cex=0.6, pos=3, offset = 0.48,  col="red") 


# Making the Figure wider, smaller text
plot(fig.matb$country.num.order, fig.matb$V2, xlab= expression(paste("Country i, ranked by ", beta[i])), ylab= expression(beta[i]))
segments(fig.matb$country.num.order, fig.matb$V3, fig.matb$country.num.order, fig.matb$V4, lwd=.5, col="gray10")
abline(h=0)
text(fig.matb$country.num.order, fig.matb$V2, fig.matb$country.nam.unorder, cex=0.5, pos=3, offset = 0.48,  col="red") 

# Checking that the labels for the countries are correct
coef(m2cyfe)
# They look correct, T&T is the lowest coefficient and it's label matches that coefficient.





### Table 6: Multilevel Model

rm(list=ls())

# This generates the same results as the old .R which used the IOworkingR.csv
data <- read.csv(file="working_2014_07_11.csv",head=TRUE,sep=",")

n <- length(data$logimports96)
y <- data$logimports96
x <- data$dispute

st1hsnum <- data$st1hs_num
country <- data$country
year <- data$year
countryyear <- data$st1_year_num

# Make moderating variable(s) by country-year
CYpolity <- data$resp_polity_it
CYlnpc <- data$lnpc
CYresp <-data$logged_resp
CYl.crisistally <-data$lcrisistally
CYcrisistally <-data$crisistally


###
#  MLM Model (in the paper): Country-Year intercepts, Country-Year Betas(dispute), 4 Mods
###
m4 <- lmer(y ~ x + CYpolity + CYlnpc + CYresp + CYl.crisistally + x:CYpolity + x:CYlnpc + x:CYresp + x:CYl.crisistally + (1 + x|countryyear) + (1|st1hsnum))
m4
#se.ranef(m4)
model4<-mtable(m4)
toLatex(model4)



